I=ones(5,1)
no_vector=151
Returns=Returns'
Innovation=Innovation'
ER=mean(Returns)*12
ECOV=cov(Returns)*12
EI=mean(Innovation)
EI=EI'
ER=ER'
INTCOV=inv(ECOV)
% compute efficient surface
a=ER'*INTCOV*ER
b=ER'*INTCOV*EI
c=I'*INTCOV*ER
d=EI'*INTCOV*EI
e=I'*INTCOV*EI
f=I'*INTCOV*I
C_determinant=a*d*f+2*b*c*e-a*e*e-b*b*f-c*c*d
d_2=(1/C_determinant)*(...
    (d*f - e*e) * INTCOV * ER ...
  + (c*e - b*f) * INTCOV * EI ...
  + (b*e - c*d) * INTCOV * I)
d_3=(1/C_determinant)*(...
    (c*e - b*f) *INTCOV * ER ...
  + (a*f - c*c) * INTCOV * EI ...
  + (b*c - a*e) *INTCOV * I)
x_0=(1/C_determinant)*(...
    (b*e - c*d) * INTCOV * ER ...
  + (b*c - a*e) * INTCOV * EI ...
  + (a*d - b*b) * INTCOV * I)
D=zeros(3,3);
D(1,1) = x_0' * ECOV * x_0;
D(2,2) = d_2' * ECOV * d_2;
D(3,3) = d_3' * ECOV * d_3;
D(1,2) = d_2' * ECOV * x_0;
D(1,3) = d_3' * ECOV * x_0;
D(2,3) = d_2' * ECOV * d_3;
D(2,1) = D(1,2);
D(3,1) = D(1,3);
D(3,2) = D(2,3);
x_1    =(1/f) * INTCOV * I
Delta_2=0.5   * INTCOV * (ER -(c/f) * I)
Delta_3=0.5   * INTCOV * (EI - (e/f) * I)
z_vector_1_1 = x_1' * ECOV * x_1;
z_vector_1_2 =ER' * x_1;
z_vector_1_3 = EI' * x_1;
figureId=1;
%<<<<<<<<<<<<<<<<<<<<< minimum-variance surface in variance <<<<<<<<<<<<<<<<<<<<<<<<<<
z_2=linspace(z_vector_1_2-2.9* z_vector_1_2, z_vector_1_2+2.9 * z_vector_1_2, no_vector);
z_3=linspace(z_vector_1_3-3.7* z_vector_1_3, z_vector_1_3+3.7 * z_vector_1_3, no_vector);
[z_2,z_3]=meshgrid(z_2,z_3);
z_1=z_2;
for i=1: no_vector,
    for j=1: no_vector,
        z_1(i,j) = [1 z_2(i,j) z_3(i,j)] * D * [1; z_2(i,j); z_3(i,j)];
    end
end
figure(figureId)
figureId=figureId+1;
surf(z_2,z_3,z_1)
axis tight
hold on
box on
%
%zlabel('Variance',  'FontName','Times New Roman', 'FontAngle','italic', 'fontsize',16.8)
zlabel('Variance',  'FontName','Times New Roman', 'FontAngle','italic', 'fontsize',12)
xlabel('Expected Return',  'FontName','Times New Roman', 'FontAngle','italic', 'fontsize',12)
ylabel('Expected GI',  'FontName','Times New Roman', 'FontAngle','italic', 'fontsize',12)
%
view([-32.5,46.0])
shading interp
%colormap(copper)
York_copper=brighten(copper,-0.23);
colormap(York_copper);
%colormap(copper);
light('Position',[ -3.19875, 0.0029, 0.97578],'Style','infinite')
%light('Position',[ -0.09875, 0.0029, 0.97578],'Style','infinite')
lighting phong
%
%MATERIAL([ka kd ks n sc]) sets the ambient/diffuse/specular strength,
%    specular exponent and specular color reflectance of the objects.
alpha(0.9)
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%
%next draw the contour (indifference curves) by contour, help contour
figure(figureId)
figureId=figureId+1;
[C,h]=contour(z_2,z_3,z_1);
set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)
hold on
% compute x
%C:\AJEJ\r_220_model_responsible_investment_by_multiple_objective\help\section_H_example\05_choose_stock__initial_compute\data_4_good_2_sin_stock_F\return__market_value__1999_2013.xlsx
%worksheet  "10 market value matrix"
%x_S=[0	0.129933007	0.137205495	0.240275264	0.492586234]';
z_S(1)=(x_S') * ECOV * x_S
z_S(2)=(x_S') * ER;
z_S(3)=(x_S') * EI;
plot3(z_S(2),z_S(3),z_S(1),'k.-','LineWidth',2,'MarkerSize',30)
xlabel('Expected return')
ylabel('Expected GI')
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%
%
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
%<<<<<<<<<<<<<<<<<<<<< nondominated surface 3D as a mesh in variance and z_S, z_SR, z_SC <<<<<<<<<<<<<<<<<<<<<<<<<<
%
%<<<<<<<<<<< next draw nondominated surface 3D as a mesh in variance <<<<<<<<<<<
no_vector=31;
z_2=linspace(z_vector_1_2-0.23 * abs(z_vector_1_2), z_vector_1_2+1.7 * abs(z_vector_1_2), no_vector);
z_3=linspace(z_vector_1_3-0.19 * abs(z_vector_1_3), z_vector_1_3+1.7 * abs(z_vector_1_3), no_vector);
[z_2,z_3]=meshgrid(z_2,z_3);
z_1=z_2;
z_4=z_1;
for i=1: no_vector,
    for j=1: no_vector,
        z_1(i,j) = [1 z_2(i,j) z_3(i,j)] * D * [1; z_2(i,j); z_3(i,j)];
        z_4(i,j) = sqrt(z_1(i,j));
    end
end
figure(figureId)
figureId=figureId+1;
mesh(z_2,z_3,z_1)
%mesh(z_2,z_3,z_4)
colormap cool
hold on
box on
alpha(0.65)
axis tight
zlabel('Variance',  'FontName','Times New Roman', 'FontAngle','italic', 'fontsize',12)
xlabel('Expected Return',  'FontName','Times New Roman', 'FontAngle','italic', 'fontsize',12)
ylabel('Expected GI',  'FontName','Times New Roman', 'FontAngle','italic', 'fontsize',12)
view([-64.0, 54.0])
%>>>>>>>>>>> next draw nondominated surface 3D as a mesh in variance >>>>>>>>>>>
%
%next compute and plot and label the criterion vector of (social) screening
plot3(z_S(2),z_S(3),z_S(1),'k.-','LineWidth',2,'MarkerSize',30)
%
%next redefine D based on
%C:\AJEJ\r_060_nondominated_surface_1x_is_1__3_objective\help\section_M_project_illustration\55_figure_40_indifference_curve\figure_40_indifference_curve.m
D(1,1) = d_2' * ECOV * d_2;
D(1,2) = d_2' * ECOV * d_3;
D(1,3) = d_2' * ECOV * x_0;
D(2,1) = D(1,2);
D(2,2) = d_3' * ECOV * d_3;
D(2,3) = d_3' * ECOV * x_0;
D(3,1) = D(1,3);
D(3,2) = D(2,3);
D(3,3) = x_0' * ECOV * x_0;
%
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
%I must manually set the following variables based on 
%C:\AJEJ\r_220_responsible_investment_by_multiple_objective\help\section_H_example\90_figure_90__N_portfolio_screen\step_a_indifference_curve__portfolio_screen
%and 
%update 20150809
%C:\new_for_this_computer\r_220_responsible_investment_by_multiple_objective\proposal\20150808_compute_end_point_N__N__zSR_zSC.doc
%
z1=z_S(1); % fixed z1 for the indifference curve
z3_for_z_H_start=10;     % higher value for locating z_H, I MUST search from higher value to lower value
z3_for_z_H_end  =0;  %  lower value for locating z_H
no_point_for_search=10000; %I need to a big number here to guarantee the precision of locating z_H and z_V
z3_step  = (z3_for_z_H_start-z3_for_z_H_end) /(no_point_for_search-1); %step for locating z_H, redefined later for indifference curve of N
%
z2_for_z_V_start=0.1031;   % higher value for locating z_V, I MUST search from higher value to lower value
z2_for_z_V_end  =0;   %  lower value for locating z_V
z2_step  = (z2_for_z_V_start-z2_for_z_V_end) /(no_point_for_search-1); %step for locating z_V
%
no_point_indifference_curve_N=500; %number of points on N
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%
%<<<<<<<<<<<<<<<<<<<<< next locate z_H <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
% the minimum-variance surface is 
% z1=[z2 z3 1] * D * ([z2 z3 1]')
% z1= D(1,1)z2^2 + D(2,2)z3^2 + D(3,3) +  2D(1,2)z2z3 +2D(1,3))z2 +2D(2,3)z3
% z1= D(1,1)z2^2  +  2 (D(1,2)z3+D(1,3)) z2  + (D(2,2)z3^2+2*D(2,3)z3+D(3,3))
% the indifference curve is
%     D(1,1)z2^2  +  2 (D(1,2)z3 +D(1,3)) z2  + (D(2,2)z3^2+2*D(2,3)z3+D(3,3)-z1) =0
%     D(1,1)z2^2  +  2 (D(1,2)z3 +D(1,3)) z2  +        (D(2,2)  z3^2 +2*D(2,3) z3 +D(3,3)-z1 )=0
% discriminant =     4*(D(1,2)*z3+D(1,3))^2 - 4*D(1,1)*(D(2,2)*(z3^2)+2*D(2,3)*z3 +D(3,3)-z1 );
%
for i=1: no_point_for_search
  z3=z3_for_z_H_start-z3_step*(i-1);
  %
  discriminant =     4*(D(1,2)*z3+D(1,3))^2 - 4*D(1,1)*(D(2,2)*(z3^2)+2*D(2,3)*z3 +D(3,3)-z1 );
  if discriminant >=0  % for the first time discriminant >=0, find z_H
    z_H(1)=z1;
    z_H(2)=                 (-2*(D(1,2)*z3+D(1,3))+sqrt(discriminant))/(2*D(1,1));
    %discard the other root (-2*(D(1,2)*z3+D(1,3))-sqrt(discriminant))/(2*D(1,1))
    z_H(3)=z3;
    %z_H
    break
  end
end
%plot3(z_H(2),z_H(3),z_H(1),'k.-','LineWidth',2,'MarkerSize',30)
%
x_H = x_0 + z_H(2)*d_2 + z_H(3)*d_3;
%next is to double-check
%x_H' * cov_matrix * x_H
%x_H' * expected_return
%x_H' * expected_CSR
%>>>>>>>>>>>>>>>>>>>>> next locate z_H >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%
%
%<<<<<<<<<<<<<<<<<<<<< next locate z_V <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
% the minimum-variance surface is 
% z1=[z2 z3 1] * D * ([z2 z3 1]')
% z1= D(1,1)z2^2 + D(2,2)z3^2 + D(3,3) +  2D(1,2)z2z3 +2D(1,3))z2 +2D(2,3)z3
% D(1,1)z2^2 + D(2,2)z3^2 + D(3,3) +  2D(1,2)z2z3 +2D(1,3))z2 +2D(2,3)z3 =z1
% the indifference curve is
% D(2,2)z3^2  +  2 (D(1,2)z2 +D(2,3)) z3  +        (D(1,1)z2^2   +2D(1,3)z2   +D(3,3)-z1) =0
% discriminant = 4*(D(1,2)*z2+D(2,3))^2 - 4*D(2,2)*(D(1,1)*(z2^2)+2*D(1,3)*z2 +D(3,3)-z1);     
for i=1: no_point_for_search
  z2=z2_for_z_V_start-z2_step*(i-1);
  %
  discriminant = 4*(D(1,2)*z2+D(2,3))^2 - 4*D(2,2)*(D(1,1)*(z2^2)+2*D(1,3)*z2 +D(3,3)-z1);
  if discriminant >=0  % for the first time discriminant >=0, find z_V
    z_V(1)=z1;
    z_V(3)=    (-2*(D(1,2)*z2+D(2,3))+sqrt(discriminant)) / (2*D(2,2));
    %discard the other root  (-2*(D(1,2)*z2+D(2,3))-sqrt(discriminant)) / (2*D(2,2))
    z_V(2)=z2;
    %z_V
    break
  end
end
%plot3(z_V(2),z_V(3),z_V(1),'k.-','LineWidth',2,'MarkerSize',30)
%
x_V = x_0 + z_V(2)*d_2 + z_V(3)*d_3;
%next is to double-check
%x_V' * cov_matrix * x_V
%x_V' * expected_return
%x_V' * expected_CSR
%>>>>>>>>>>>>>>>>>>>>> next locate z_V >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%
%
%<<<<<<<<<<<<<<<<<<<<< next locate z_SR <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
if z_S(3)>z_V(3) 
  %zS is above the broken line (z3=z3 of zV),
  z_SR(1)= z_S(1);
  z_SR(3)= z_S(3);
  %     D(1,1)z2^2  +       2 (D(1,2)z3 +D(1,3)) z2       +        (D(2,2)  z3^2      +2*D(2,3) z3      +D(3,3)-z1 )=0
  % discriminant =          4*(D(1,2)*z3+D(1,3))^2      - 4*D(1,1)*(D(2,2)*(z3^2)     +2*D(2,3)*z3      +D(3,3)-z1 );
    discriminant =          4*(D(1,2)*z_SR(3)+D(1,3))^2 - 4*D(1,1)*(D(2,2)*(z_SR(3)^2)+2*D(2,3)*z_SR(3) +D(3,3)-z_SR(1));
  z_SR(2)=                (-2*(D(1,2)*z_SR(3)+D(1,3))+sqrt(discriminant)) / (2*D(1,1));
  %discard the other root (-2*(D(1,2)*z_SR(3)+D(1,3))-sqrt(discriminant)) / (2*D(1,1))
else
  z_SR=z_V;
end
%next plot and label z_SR
plot3(z_SR(2),z_SR(3),z_SR(1),'k.-','LineWidth',2,'MarkerSize',30)
%next compute x_SR (the weight), based on x = x0 + e2d2 + e3d3, model (21) paper_r_060.pdf
%C:\AJEJ\r_060_nondominated_surface_1x_is_1__3_objective\version_20xxxxxx\20141024_finish_fourth_draft_Steuer_appraisal
x_SR = x_0 + z_SR(2)*d_2 + z_SR(3)*d_3
z_SR;
%next is to double-check
%x_SR' * cov_matrix * x_SR
%x_SR' * expected_return
%x_SR' * expected_CSR
%>>>>>>>>>>>>>>>>>>>>> next locate z_SR >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%
%
%<<<<<<<<<<<<<<<<<<<<< next locate z_SC <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
if z_S(2)>z_H(2) 
  %zS is to the right of the broken line (z2=z2 of zH)
  z_SC(1)= z_S(1);
  z_SC(2)= z_S(2);
  %     D(2,2)z3^2  +       2 (D(1,2)z2      +D(2,3)) z3  +          (D(1,1)z2^2        +2D(1,3)z2       +D(3,3)-z1) =0    
  % discriminant =          4*(D(1,2)*z2     +D(2,3))^2   - 4*D(2,2)*(D(1,1)*(z2^2)     +2*D(1,3)*z2     +D(3,3)-z1)
    discriminant =          4*(D(1,2)*z_SC(2)+D(2,3))^2   - 4*D(2,2)*(D(1,1)*(z_SC(2)^2)+2*D(1,3)*z_SC(2)+D(3,3)-z_SC(1));
  z_SC(3)=                (-2*(D(1,2)*z_SC(2)+D(2,3))+sqrt(discriminant)) / (2*D(2,2));
  %discard the other root (-2*(D(1,2)*z_SC(2)+D(2,3))-sqrt(discriminant)) / (2*D(2,2))
else
  z_SC=z_H;
end
%next plot and label z_SC
plot3(z_SC(2),z_SC(3),z_SC(1),'k.-','LineWidth',2,'MarkerSize',30)
%next compute x_SC (the weight), based on x = x0 + e2d2 + e3d3, model (21) paper_r_060.pdf
%C:\AJEJ\r_060_nondominated_surface_1x_is_1__3_objective\version_20xxxxxx\20141024_finish_fourth_draft_Steuer_appraisal
x_SC = x_0 + z_SC(2)*d_2 + z_SC(3)*d_3
z_SC;
%next is to double-check
%x_SC' * cov_matrix * x_SC
%x_SC' * expected_return
%x_SC' * expected_CSR
%>>>>>>>>>>>>>>>>>>>>> next locate z_SC >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%
%
%<<<<<<<<<<<<<<<<<<<<< next locate indifference curve of N <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
% z3 varies from z_V(3) to z_H(3)
z3_step  = (z_H(3)-z_V(3)) /(no_point_indifference_curve_N-1);
k=0;
for i=1: no_point_indifference_curve_N
  z3=z_V(3)+z3_step*(i-1);
  discriminant =                4*(D(1,2)*z3+D(1,3))^2 - 4*D(1,1)*(D(2,2)*(z3^2)+2*D(2,3)*z3 +D(3,3)-z1 );
  if discriminant >=0  % discriminant >=0, find a point
    k=k+1;
    indifference_curve_N(k,1)=z1;
    indifference_curve_N(k,2)=(-2*(D(1,2)*z3+D(1,3))+sqrt(discriminant))/(2*D(1,1));
    %discard the other root   (-2*(D(1,2)*z3+D(1,3))-sqrt(discriminant))/(2*D(1,1))
    indifference_curve_N(k,3)=z3;
  end
end
%
%next plot the indifference curve
plot3(indifference_curve_N(:,2),indifference_curve_N(:,3),indifference_curve_N(:,1),'k','LineWidth',2)
%>>>>>>>>>>>>>>>>>>>>> next locate indifference curve of N >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%
%
%>>>>>>>>>>>>>>>>>>>>> nondominated surface 3D as a mesh in variance and z_S, z_SR, z_SC >>>>>>>>>>>>>>>>>>>>>>>>>>
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% compute out-of-sample
SC_Out_Return=Out_Return*x_SC
SC_Out_GI=Out_GI*x_SC
SR_Out_Return=Out_Return*x_SR
SR_Out_GI=Out_GI*x_SR
VW_Out_Return=Out_Return*x_S
VW_Out_GI=Out_GI*x_S
Return_vector=table(SC_Out_Return,SR_Out_Return,VW_Out_Return)
GI_vector=table(SC_Out_GI,SR_Out_GI,VW_Out_GI)
% compute portfolios in-between
SF_Out_Return=Out_Return*x_SF
SF_Out_GI=Out_GI*x_SF
